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Abstract. We present a numerical approach for solving the free bound- 
ary problem for the Black- Scholes equation for pricing American style 
of floating strike Asian options. A fixed domain transformation of the 
free boundary problem into a parabolic equation defined on a fixed spa- 
tial domain is performed. As a result a nonlinear time-dependent term 
is involved in the resulting equation. Two new numerical algorithms are 
proposed. In the first algorithm a predictor-corrector scheme is used. The 
second one is based on the Newton method. Computational experiments, 
confirming the accuracy of the algorithms are presented and discussed. 



1 Introduction 

In this paper we consider the problem of pricing American style Asian options, 
analyzed by Bokes and the second author in [JJ (see also [JTJ). Asian options be- 
long to the group of the so-called path-dependent options. Their pay-off diagrams 
depend on the spot value of the underlying asset during the whole or some part(s) 
of the life span of the option. Among path-dependent options, Asian option de- 
pend is on the arithmetic or geometric average of spot prices of the underlying 
asset. During the last decade, the problem of solving the American option prob- 
lem numerically has been subject for intensive research [11619110(13] (see also 
[TT] for overview). A comprehensive introduction to this topic can be found in 
[BJ. Comparison of various analytical and numerical approximation methods of 
calculation of the early exercise boundary a position of the American put option 
paying zero dividends is given in [7. . An improvement of Han and Wu's algorithm 
[4] is described in [14] . Our goal is to propose and investigate two front-fixing 
numerical algorithms for solving free boundary value problems. The front-fixing 
method has been successfully applied to a wide range of applied problems arising 
from physics and engineering, cf. [318] and references therein. The basic idea is 
to remove the moving boundary by a transformation of the involved variables. 
Transformation techniques were used in the analysis and numerical computa- 
tion of the early exercise boundary in the context of American style of vanilla 
options [JO] as well as Asian floating strike options [1111112] , In comparison to 



the existing computational method [I] we do not replace the algebraic constraint 
by its equivalent integral form (see [1112] for details) which is computationally 
more involved. In this paper we solve the corresponding parabolic equation with 
an algebraic constraint directly as it was proposed in The approach pre- 
sented in however suffered from the necessity of taking very small time 
discretization steps. Here we overcome this difficulty by proposing two new nu- 
merical approximation algorithms (see Section 4). They are based on the novel 
technique proposed by the first author and Valkov in [5]. We extend this ap- 
proach for American style of Asian options. In Section 5, a numerical example 
illustrating the capability of our algorithms are discussed. 

2 The Free Boundary Problem 

Following the classical Black-Scholes theory, the second author and Bokes [I] 
analyzed the problem of pricing Asian options with arithmetically averaged strike 
price by means of a solution to a parabolic PDE with a free boundary Sf(t, A): 
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<t <T, < S < Sf(t,A), satisfying the boundary conditions 

V(t, 0, A) = 0, for any A > and < t < T, (2) 
^(t,S f (t > A),A) = l, V(t,S f (t,A),A) = S f (t,A)-A, (3) 

and the terminal condition (terminal pay-off condition) at the maturity time T: 

V(T, S, A) = max(5 - A, 0), S,A>0. (4) 

Here S > is the stock price, A > is the averaged strike price, r > is 
the risk-free interest rate, q > is a continuous dividend rate and a > 
is the volatility of the underlying asset returns. The arithmetically averaged 
price A = A t calculated from the price path {S u ,u e [0, T]} at the time 
t is defined as At = \ L S u du. For floating strike Asian options, it is well 
known (see e.g. [6 2 I ) that one can perform a dimension reduction by intro- 
ducing a new time variable r = T — t and a similarity variable x defined as: 
x = A/S, W(x, t) = V(t, 5, A) J A. The spatial domain for the reduced equa- 
tion is given by l/p(r) < x < oo, r e (0,T), p(r) = Sf(T - t,A)/A. Following 
([10 13 1), we can apply the Landau fixed domain transformation for the free 
boundary problem by introducing a new state variable £ and an auxiliary func- 
tion iT(^,r) = W(x, r) + x^-(x, r), representing a synthetic portfolio. Here 
£ = In (p(r)x). In [1110113] it is shown that under suitable regularity assump- 
tions on the input data the free boundary problem (|TJ)— ((4j) can be transformed 



into the initial boundary value problem for parabolic PDE: 
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The coefficients a and ft are defined as follows: 

a(?>TJ = -r^ + r-g- , /3(£,r)=r + - . (7 

p(Tj 2 1 — T 1 — T 

According to pQ the free boundary function p(r) and the solution 77 should fulfill 
the constraint: 

P(r) = / ,V / » P(0) = max -L— ,1 . (8) 



l + q{T-T) ' ^ v ' \l + qT 

As for derivation of the initial free boundary position p(0) in © we refer to 
PQ or [50. A solution 77 to the problem ©-© is continuous for t > 0. The 
discontinuity appears only at the point P* = (ln(/?(0)), 0). The derivatives of 
the solution exist and are sufficiently smooth in [0, L] x [0, T), outside of a small 
neighbourhood of P*. Another important fact to emphasize is that for times 
t — > + (i.e. when r — > T) the coefficients a, ft become unbounded. 



3 Finite Difference Schemes 

In order to solve the problem ©-© numerically, we introduce L which is suffi- 
ciently large upper limit of values of the £ variable (a safe choice is to take L is 
equal to five times ln(p(0))), where we prescribe 77(L, r) = 0. Next, for given pos- 
itive integers N and M we define the uniform meshes: ZJh = {0}U{L}UuJh, u>h = 
{& =ih, i = 1, . . . , (TV — 1), h = L/N} and ZJ k = {0} U {T} U u k , uj k = fa = 
jk, j — 1,...,(M — 1), k — T/M}. Our goal is to define a finite difference 
method which is suitable for computing y\ m II (£i,7j-) for (£i,r,-) G uih x <^k and 
associated front position z 1 sa pfa) for Tj G ujk- The implicit difference scheme 
has the following form: 

nP V j+1 - V j+1 rr 2 V j+1 - 2v j+1 + lP + 1 

2/i 2 h 2 +P 2/4 ~ A ' 



(10) 
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For the initial condition for the free boundary we have z° = p(0). An algebraic 
nonlinear system of equations can be derived from ([3]) for i = 1, . . . , N — 1 , (|10l) 
and (|12p. In [S] the authors apply implicit finite difference scheme, semi-implicit 
scheme and upwind explicit scheme for the American put option, combining with 
the penalty method. The time step parameter for the explicit case is very small, 
k = 5.0 ■ 10~ 6 . Therefore in this work we consider the case of a fully implicit 
scheme. One can also apply a scheme of the Crank-Nicolson type. 



4 Numerical Algorithms 

In order to solve the nonlinear system of algebraic equations we developed the 
following two algorithms. 

Algorithm 1. This algorithm is based on the predictor- corrector scheme and 
consists in the following steps, (see also [15116] for the case of pricing American 
put options). 

Step 1. Predictor. Let the solution and the free boundary position on the 
time level tj be known. Instead of (|12p we use another approximation of (J5J by 
introducing an artificial spatial node 

(1 + q(T r J+1 )) zJ +1 = 1 + r(T - r j+1 ) + y (T - r J+1 ) V{+ ~^ . (13) 

An additional equation can be obtained from ([5]) by taking the limit £ — > and 
using the fact that d T II(0, t) — 0: 

*i +l Vl ^ Y~ % +/P + Y +1 - °- ( 14 ) 

Using (fT3|) we can express y^ 1 as: 

ylX 1 = y{ +1 ( - r + % (15) 

1 -T j+X J ®~ 



Inserting it into (|14p we conclude the following equation for the value y{ +1 : 



«r = + 7 ) - r +¥ ^ 7 J - - 1. (16) 

Instead of the implicit scheme © we make use of its explicit variant for i = 1 
in order to derive 

yJ^_A + y\_A _ £ d^M + ri + = Q (17) 

This way we obtain a nonlinear system (p~6|) . (fTT|) for unknowns and z J+1 . 
The system is indeed nonlinear as a\ +1 depend on Now, by replacing 

y{ +1 O and O P +1 we construct the predictor value of P +1 . 



if 



Step 2. Corrector. We again use Equation @ in a slightly different form: 
i+i 



+ a 
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P +1 yi +1 =0, (18) 



where approximation a{ +1 takes into account the already constructed predictor 
value i.e. 



r-Q~ T 



P' +1 exp(-&)-l 



T-t. 



3+1 



(19) 



Next we use the corrected solution y{ +1 and Equation (IT^I) in order to obtain 
the corrected value for the free boundary position on the next time layer. 
Algorithm 2. We now describe an algorithm based on the Newton method. A 
variant of this method was applied for an American Call option problem in [5]. 

Step 1. We eliminate the known boundary values j/q +1 = —1 and y^ 1 = 
from J9J. Taking into account (fT2|) we obtain a nonlinear system for N unknowns: 



J+ 1 



i 



y\' , t = 1,2, N— 1 and We denote by Y the vector of these N unknowns 
at the l-th iteration. 



/ i I i i 

Step 2. We have to solve the equation F= with F= ( Fi F2 



where 



Fiji = 1)2, correspond to Equations ([9]) and (|12|) . respectively. To this end, we 

apply the Newton method in the following form: J ( Y — Y)= — F, with the 

1 1 

Jacobi matrix defined by: J= {Jij)i,j=i,2 where 
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J11 
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a J+1 b J+1 

u N-2 c N-2 u N-2 
J+ 1 J+ 1 
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9a W-2„,j + l 



0b W + -2„J + l 



a 2 j+i j/jv-2 
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J 2i=(7^r ; 4k' 0,...,0) where D = q + 1/(T -t-> +1 ) and J 22 = 1. Similarly 

T 



^YiY 2 ^ , Yi= ...,yj^\y Y 2 = As for the elements of the 



matrix Jn we have: 



and d 3 { +1 = l/(2h)(z j+1 exp(— £j) — 1)/(T — T j+i)- The iteration process is re- 

i+i i 

peated until the condition || Y — Y || < tol is fulfilled. 

Step 3. The solution on the (j + 1 )-th time layer is considered as an initial 

i i+i 1 i 

iteration for the next time layer. For solving J ( Y — Y)= — F we perform the 

i i+i i 

following stages. First, we solve the linear system of equations Jn Y i= — Fi 

i i i i+i i i i 
+ JnYi — J12 Y 2 + Ji2Y2- Since the matrix Jn is tridiagonal we can apply 

i+i i i+i i i+i i 

the Thomas algorithm to find Y i- Next, we solve J i2 Y i + J22 Y 2= — F2 ■ 

Remark 1. In both algorithms we choose the last time step k — e with e = 10~ 7 , 
i.e. tm = T — e. To overcome possible numerical instabilities of these methods for 
r — > T (i.e. t — ?- 0) we use the so called upwind and downwind approximations of 
the term - — j^ P 7 f.~^' ) ~ 1 ^ depending of the sign of the term exp(— £j) — f . 

5 Numerical Experiments 

In this section we consider problem ([T]) with parameter values r = 0.06, q = 0.04, 
a = 0.2 and T = 50, taken from examples presented in pQ. Since there exists 
no analytical solution to the proposed free boundary problem, we use the mesh 
refinement analysis with doubling the mesh size h. In Tab. 1 we present the 
position of the free boundary position p(r) at different times r constructed by 
the Newton method. We also present the difference between two consecutive 
values and the convergence ratio are presented. The results show nearly first 
order of accuracy for the free boundary and the CR increases with increasing 
r (see Tab. 1). In Fig. [Qi) a 3D plot of the portfolio function 77 for T = 50, 
N = 200, M = 500 is presented. In Fig. [Qd) the profiles of the function 77(£ , r) 
for t = 0, 0.1, 10, 25, 50 obtained by the Newton method are depicted. 

In Fig. 2a) we show a comparison of the free boundary position p(r) computed 
by our two algorithms (Predictor-corrector and Newton's based method) and by 
numerical methods from [I] (Bokes) and [2] (Kwok) . It turns out that the New- 
ton's based method gives nearly the same results as those of [112] . On the other 
hand, predictor-corrector methods slightly underestimates the free boundary po- 
sition p(r). In Fig. 2b) we show the free boundary position Xf(t) = l/p(T — t) 
for the original model variables x = A/S and t. The continuation region and 
exercise region are also indicated. 

6 Conclusions 

In this paper we have analyzed numerical algorithms for solving the free bound- 
ary value problem for American style of floating strike Asian options. To solve 
corresponding degenerate parabolic problem we have applied Landau's front fix- 
ing transformation method. We proposed two numerical algorithms based on the 
predictor-corrector scheme and the Newton's method. The predictor-corrector 



Table 1. Mesh-refinement analysis and the convergence ratio (CR) of the New- 
ton method. 



N 


P (t = 10) 


difference 


CR 


p(r = 20) 


difference 


CR 


p(r = 40) 


difference 


CR 


50 


1.949988 






1.991675 






1.796663 






100 


1.955552 


5.5640c-3 




1.995525 


3.8502c-3 




1.803276 


6.6133e-3 




200 


1.958037 


2.4850c-3 


1.16 


1.996945 


1.4194e-3 


1.44 


1.805149 


1.8729e-3 


1.82 


400 


1,959199 


1.1617o-3 


1.10 


1.997515 


5.7099o-4 


1.31 


1.805667 


5.1799c-4 


1.85 


800 


1.959758 


5.5965e-4 


1.05 


1.997765 


2.4919o-4 


1.20 


1.805813 


1.4621e-4 


1.82 




a) b) 

Fig. 1. (a) A 3D plot of the portfolio function 77 for T = 50, N = 200, M = 500; 
(b) Profiles of the function I7(f , r) for r = 0, r = 0.1, r = 10, r = 25, t = 50. 

scheme is computationally faster when compared to the Newton method. It yields 
a good approximation close to expiry. However, its accuracy is decreased for 
times close to the initial time. The second algorithm based on Newton's method 
yields better approximation results over the whole time interval. Although all 
finite difference approximations are of second order, due to discontinuity of the 
initial datum and nonlinear behavior of the coefficients in all discrete equations, 
the results show nearly the first order rate of convergence. 
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